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CN Abstract 

o 

£f) Different numerical approaches for the stray-field calculation in the context of micromag- 

netic simulations are investigated. We compare finite difference based fast Fourier trans- 

■^j- form methods, tensor grid methods and the finite- element method with shell transformation 

in terms of computational complexity, storage requirements and accuracy tested on several 
benchmark problems. These methods can be subdivided into integral methods (fast Fourier 
transform methods, tensor-grid method) which solve the stray field directly and in differen- 
J> tial equation methods (finite-element method) , which compute the stray field as the solution 

ITj of a partial differential equation. It turns out that for cuboid structures the integral meth- 

ods, which work on cuboid grids (fast Fourier transform methods and tensor grid methods) 
outperform the finite-element method in terms of the ratio of computational effort to accu- 
racy. Among these three methods the tensor grid method is the fastest. However, the use 
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of the tensor grid method in the context of full micromagnetic codes is not well investigated 
yet. The finite-element method performs best for computations on curved structures. 

Keywords: micromagnetics, stray-field, fast Fourier transform, tensor grid methods, low- 
rank magnetization, finite- element method 

1 Introduction 

Micromagnetic simulations nowadays are highly important for the investigation of ferromagnetic 
materials which are used in storage systems and electric motors and generators. In these 
simulations the magnetic state of the ferromagnet is represented by a classical magnetization 
vector field. 

The computation of the non-local magnetostatic interactions is thereby the most time- 
consuming part. Naive implementation of the superposition-based integral operators (5) or 
solvers for the underlying differential equation (Poisson equation (3)) yield computational costs 
proportional to the square of the number of grid points, i.e. 0{N 2 ). Several methods have been 
introduced in literature in order to reduce these costs. 

The magnetic scalar potential can be computed by solving the Poisson equation. The solu- 
tion of the Poisson equation with the finite-element method (FEM) has a complexity of 0{N) if 
boundary conditions are given at the boundary of the sample and a multigrid preconditioner is 
used [28]. However, the stray-field problem has open boundary conditions, where the potential 
is known at infinity. Two possible solutions for the open boundary problem are the coupling of 
the boundary element method (BEM) with the finite-element method [16] and the application of 
a shell transformation [6]. BEM gives an additional complexity of 0(M 2 ) where M is the num- 
ber of boundary nodes. This complexity can be reduced to O(MlogM) by application of the 
^-matrix approximation for the dense and unstructured boundary element matrices [7, 21, 27]. 
The storage requirements and computational complexity of the FEM with shell transformation 
will be described in the forthcoming text. 

Another class of methods rely on the evaluation of volume and/or surface integrals for the 
direct computation of the magnetostatic potential or the field, e.g. fast multipole methods [4, 5], 
nonuniform grid methods [23] and fast Fourier transform (FFT) methods [2, 25], scaling from 
0{N) to 0(N\ogN). The more recent tensor grid method (TG), which also belongs to this 
class scales even better under certain assumptions. 

In this paper we compare recently developed algorithms, namely the FFT-based methods for 
the computation of the field via the scalar potential (SP) and directly (DM) [1, 2, 11], a recently 
developed approach from numerical tensor-structered methods (TG) [13], and the finite-element 
method with shell transformation (FES), which is a FEM method that does not rely on BEM 
approaches and thus only introduces sparse matrices. 



Consider a magnetization configuration M that is defined on a finite region £1 = {r : M(r) ^ 
0}. In order to perform minimization of the full micromagnetic energy functional or solve the 
Landau-Lifshitz-Gilbert (LLG) equation it is necessary to compute the stray field within the 
finite region f2. The stray- field energy is given by 



2 Stray-Field Problem 




(i) 
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The Landau-Lifshitz-Gilbert equation reads 

M t = ~-^M x H cS + "1 2 , M x (M x H eS ), (2) 
1 + a z M s (l + or) 

where a is the Gilbert damping constant and H e Q is the effective field given by the negative 
variational derivative of the energy density. In both cases the stray field is only required to be 
known within Q. The stray field H has a scalar potential <j>, which is the solution of a Poisson 
equation [19] 

H = -V<f> (3) 
Acp = V • M (4) 

The stray field H and thus also the scalar potential (f> are required to vanish at infinity. This 
boundary condition is often referred to as open boundary condition [15]. 



3 Methods 

3.1 FFT Methods (SP and DM) 

One way to reduce the computational complexity is to employ the fast Fourier transform (FFT). 
FFT methods solve an integral solution of the Poisson equation by applying the convolution 
theorem. The solution to the Poisson equation (3) is given by the integral, see [19], 

1 f V ■ M(r') j3 , 1 f n'-M(r') JAl 
4vr7 n \r-r'\ 4vr J m \r - r'\ 

which directly fulfills the required open boundary condition. Performing integration by parts 
yields 

Mr) = -L f M(r') • V / 1 — ^d 3 r' (6) 

= S(r -r')*M(r'). (7) 
By employing the convolution theorem 

cj) = S*M = F- 1 (f{S)-F{M)], (8) 



this convolution can be discretized and solved with the fast Fourier transform. A prerequisite for 
this procedure is the usage of an equidistant grid, which is required for a discrete convolution. 
The stray field 

H(r) = -Vcf>(r), (9) 

can be obtained by applying finite differences. This method is referred to as the scalar-potential 
method (SP) in the following. It is described in detail in [2]. 

It is also possible to compute the stray field H directly as a result of a matrix-vector 
convolution. 

H(r) = -~ [ (VV . 1 W(r')dV (10) 



Air J \ \r 

= N{r-r')*M{r'). (11) 

Here N denotes the demagnetization tensor. Similar to (8) the convolution can be solved as an 
element-wise matrix-vector multiplication in Fourier space. This method is referred to as the 
demagnetization-tensor method (DM) in the following and is implemented by different finite- 
difference codes [10, 11, 29]. For the numerical experiments we use MicroMagnum [1, 11] which 
implements both the SP and the DM method. 
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grid spacing 

< > 

midpoint spacing 
Figure 1: Grid spacing and midpoint spacing in TG Methods 

3.2 Tensor Grid Methods (TG) 

Tensor grid methods (TG) for micromagnetic stray-field computation were recently introduced 
in [13] and [17]. They were developed for the purpose of handling so called low-rank tensor or 
compressed tensor magnetization, see [22] for a survey, in order to accelerate the computations 
and relieve storage requirements, see [12]. In the following we give a brief introduction into the 
ideas behind this method, also see [13] for a detailed description. 



3.2.1 Analytical preparations 

The computation of the stray field within the magnetic body is based on the explicit integral 
formula for the scalar potential (6). The main idea is the usage of a representation for the 
integral kernel in (6) as an integral of a Gaussian function by the formula 

2 /^e-Vl— 'P dT> (12) 
\r-r'\ 6 vWr 

which leads from (6) to 

</>( r ) = _L f T 2 f e - T ^ r - r '^M(r') ■ (r -r')d 3 r' dr. (13) 
27T2 Jr Jn 

Equation (13) reduces the computation to independent spatial integrals along each principal 
direction (the part of the O integral without the magnetization is now a product of independent 
1-D integrals). This analytical preparation directly results in a reduction of the computa- 
tional effort from 0(N 2 ) to C(iV 4 / 3 ) if discretized on a tensor product grid before even using 
compressed/low-rank tensor formats for the discretized magnetization components. A similar 
method was introduced for the computation of the electrostatic scalar potential, [20]. 
The additional r-integration is carried out by the exponentially convergent Sine quadrature [18], 
the spatial integrals are computed by Gauss-Legendre quadrature, both resulting in a numerical 
error of about the machine epsilon. 



3.2.2 Discretization on a tensor product grid 

The magnetic body is discretized on a tensor product grid arising from the tensor outer 
product of three vectors h p G R p ,p = 1 ... 3 related to the grid spacings along each axis (see 
Fig. 1). This results in a not necessarily uniform Cartesian grid but in contrast to methods like 
DM/SP described before, tensor grid methods make use of the tensor-product interpretation of 
such grids. 

The magnetization on the center points of the cells is given as a 3-tensor [22] for each component, 
i.e. 

€R Nl * NaxNa , p=1...3 (14) 
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where Ni, N2, N3 are the number of cells in the principal directions. Thereby it is possible to use 
low-rank representation for the magnetization like Canonical/Parallel Factors Decomposition 
(CP) or Tucker formats, see the Appendix A or [22]. We obtain the potential on the center 
points of the computational cells, as the discrete analogue of (13), by a so-called block- CP tensor 
[13] 



3 R 

vNixN 2 xN 3 

■>-■'< - 

P=l 1=1 



1 3 R 

3 * = 2^ ^5> ; sinh(r z ) 2 M^ x 1 D[x 2 D^x 3 D l z . (15) 



Gaussian matrices Dl G R n i xN i come from 



Here (77, w/) are the nodes and weights arising from the Sinc-quadrature with R terms. The 

come from 

d L,j q := I g(xt,x',Ti)dx', (16) 



"i 9 



D\--={4 qk )- (17) 

where Qj q denotes the j-th interval on the (partitioned) g-th axis with length (h q )j and x\ is 
the midpoint of the i-th interval on the g-th axis. The Gaussian integrands are given by 



■ exp(-sinh(r) 2 (a-a') 2 ) QJ^P, 
9 (g> (a,a,T) :={ . , . 2 , 2 (18) 

(a - a J exp(— sinn(r) (a — a ) J g = p, 



and are approximated using Gauss-Legendre quadrature, as mentioned above. The field within 
0, is derived from (9) by finite-difference operators of second order. 

3.2.3 Low-Rank Magnetization 

Equation (15) allows the treatment of specially structured magnetization tensors, like CP or 
Tucker tensors [22] that have a reduced number of degrees of freedom, see Tab. 1, and accelerate 
the computation up to sub-linear effort (below the volume size iV 3 ), see Tab. 2. As a conse- 
quence TG methods using low-rank magnetization allow larger models with finer discretization 
density than conventional methods. 

We now show by means of numerical experiments that typical single domain states [26] have 
highly accurate low-rank representations. Fig. 2 shows the approximation properties of a flower 
and a vortex state as described in Sec. 5 via the CP format and the Tucker format using an 
alternating least squares algorithm (ALS) [22] for the approximations. The plots 2a and 2b are 
computed independently from random initial guesses used in the ALS algorithm. We set the 
parameters in (27) as a = c = 0.5, 6=1 and in (28) as r c = 1/2. Fig. 2a shows the dependence 
of the relative error (19) w.r.t. the rank for fixed discretization density, where Fig. 2b indicates 
the dependence w.r.t. the discretization density iV for fixed rank. 
The relative errors are measured in the Frobenius norm, i.e. 



dense "^^low- 
p=x,y,z p=x,y,z 



relerr=( £ llM^ - M£. rank " ) a /( £ \\^L IV ■ (19) 



The Tucker format generally leads to a better approximation, where Fig. 2b essentially shows 
no loss of accuracy while increasing the mesh density. 
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Figure 2: Low-rank approximation of flower and vortex state via Tucker and CP decomposi- 
tion, (a) Relative error w.r.t. approximation rank r. N = le+06. (b) Relative error w.r.t. 
discretization density N. Rank r = 5. 



(a) 



(b) 





Figure 3: Parallelepipedic shell surrounding the cuboid sample, (a) The transformation is 
carried out along the blue line. The origin of the transformation moves along the yellow middle 
plane, (b) Since the area of interest is the sample, the mesh is coarsened in the shell. 



3.3 FEM Methods (FES) 

Within the finite-element framework the Poisson equation (3) is solved by the weak formulation 



V<f> ■ Vv <Tx 



M-Vvd d x V v,(j)£V 



(20) 



where Dirichlet boundary conditions are embedded in the trial function space V, i.e. the function 
space of the solution (p. In the case of the stray-field problem the boundary conditions at the 
sample boundary dVL are unknown. They are defined as zero at infinity. The treatment of these 
open boundary conditions is the main difficulty for finite-element stray-field calculations. 

We present the results for a transformation technique. The sample is surrounded by a finite 
shell which is also meshed. A bijective transformation from the finite shell to the complete 
exteriour of the sample is applied by introducing a metric tensor to the weak formulation. The 
particular transformation we use is known as "parallelepipedic shell transformation" [6]. The 
sample is put into a cuboid volume and a shell consisting of six parallelepipeds is created (see 
Fig. 3a). 
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Figure 4: Sketch of the shell transformation in two dimensions. Within each shell patch the 
shell points are transformed in a radial sense. In order to achieve continuity between the patches 
the origin O of the one-dimensional transformation has to be continuous between the patches. 
We choose the origin O to move on the middle plane of the sample. The third dimension is 
treated in the same manner. 



The transformation is chosen such that points located at the inner boundary of the shell are 
mapped to themselves. Points on the outer boundary of the shell are mapped to infinity. The 
Jacobian of the transformation is requested to be 1 on the inner boundary of the shell in order 
to be continuous across the sample boundary. 

These conditions still leave some space for the choice of transformation. The most important 
aspect of this method is the distortion of the test and trial functions in the transformed area. 
In order to get a good result, the test and trial functions must be distorted such that they 
are able to model the natural decay of the magnetic potential. This obviously also depends on 
the choice of test and trial functions. From (6) it is seen that the potential decays with 1/r 2 . 
We choose our test and trial function 4>h-, v h G Vh to be continuous and piecewise third-order 
polynomial (V-s) 

Vh = {v h G H\0) : v h \ T G V 3 (T) V T G T h ) (21) 

where H 1 is a Sobolev space and Th is a tetrehedron tesselation (see Fig. 3b). The transformation 
per shell patch is carried out in a radial sense as sketched in Fig. 4. The scalar transformation 
is given by 

* = *'S < 22 » 

with Rx, i?2, x, and X as shown in Fig. 4. This transforms the third order polynomial test and 
trial functions as 

a + bx + cx 2 + dx 3 ->• a' + bf— + c'-^o + d '4v ( 23 ) 

X X X 

The discretized weak formulation then reads 

[ (VcPhfgVvh d 3 x = ( M ■ Vv h d 3 x V v h eV h (24) 
Jn Jn 

{1 if X £ ^sample 

j-lT . . i (25) 

J I J I J II X G iishell 

where r2 samp i e and figheii denote the disjoint regions of the sample and the transformed shell 
with ^sample U ^sheii = ^ an d J is the Jacobian matrix of the transformation. This directly 
translates to a linear system of equations, where the solution vector contains the coefficients in 
terms of the discrete function basis. The size of this solution vector is referred to as degrees of 
freedom (DoF). The implementation of this method is done with FEniCS [24]. 
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Method Setup Magnetization Potential Field 
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ps 48iV 


3N 


AT 
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Table 1: Storage in number of floating point values w.r.t. number of computational cells/degrees 
of freedom N. In TG methods r denotes the tensor rank and R denotes the number of Sinc- 
quadrature nodes. 

4 Storage Requirements and Computational Complexity 

The costs of the different methods are compared in terms of storage requirements and com- 
putational complexity. Table 1 and 2 show the results. We choose iV to be the number of 
computational cells in the case of DM, SP and TG methods. In the case of finite-element 
methods (FES) N refers to the number of degrees of freedom. 

Besides the memory needed for the storage of the magnetization configuration and the 
stray field, all methods require a certain amount of extra storage for auxiliary constants. In 
case of the DM and SP methods this includes the convolutions kernels, TG needs the one 
dimensional Gaussian matrices, and finite-element methods (FES) require the stiffness matrix 
as an auxiliary constant. These constants depend on the geometry and discretization only. 
This means that the computation of auxiliary constants has to be done only once for different 
magnetization configurations. Thus their complexity is almost irrelevant in the context of LLG 
computations and energy minimization. The storage requirements for these constants as well as 
the computational complexity of their calculation is summarized in the setup column in both 
tables. 

Storage requirements for TG methods depend on the rank r used for the low-rank tensor 
representation of the magnetization components. Often the rank is much smaller than iV 1 / 3 , 
the discretization size in one spatial dimension. This makes the storage requirements for the 
magnetization, potential and field proportional to rN 1 / 3 . For the setup the (iV 1 / 3 x A" 1 / 3 ) 
Gaussian matrices need to be computed and stored, thereby R in Tab. 1 and 2 denotes the 
number of Sinc-quadrature nodes. The computational effort in TG methods also depends on the 
tensor format used for the representation of the magnetization, see Tab. 2. If the magnetization 
has a low-rank representation, TG methods usually reduce this complexity below the number 
of computational cells (sub-linear), making this methods the fastest available nowadays. 

The storage requirements for the other three methods are proportional to N, which is a 
result of the dense representation of the magnetization, see Tab. 1. A well-known result is the 
NlogN complexity of the convolution in FFT methods (DM/SP), likewise this is the asymptotic 
operation count for those methods, Tab. 2. In the FES method sparse linear systems have to 
be solved for the computation of the scalar potential. We used a conjugate gradient solver 
(CG) with an algebraic multigrid preconditioner (AMG) and measured the complexity w.r.t. A" 
(DoF) experimentally, finding a linear dependence on the system size (with a small logarithmic 
scaling factor). 
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Method Setup Potential Field Energy 
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Table 2: Computational complexity w.r.t. number of computational cells/degrees of freedom 
N . In TG methods r denotes the tensor rank. Every column depends on its left neighbor, e.g. 
the calculation of the field requires the previous calculation of the potential etc. 



(a) (b) (c) 




Figure 5: Magnetization configurations in a 1 x 1 x 1 cube used for numerical experiments. The 
magnetization is normalized, its direction is color coded, (a) homogeneous magnetization (b) 
fower state (c) vortex state. 



(a) (b) 




Figure 6: Magnetization configurations in a 1 x 1 x 0.1 cuboid used for numerical experiments. 
The magnetization is normalized, its direction is color coded, (a) homogeneous magnetization 
(b) vortex state. 
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Figure 7: Convergence of the calculated stray-field energy for different geometries and magneti- 
zation configurations. Like in Tab. 1 and 2, N is the number of cells for the tensor grid methods 
(DM, SP and TG) and the number of degrees of freedom in the case of finite elements (FES), 
(a) homogeneously magnetized lxlxl cube (b) flower state in 1 x 1 x 1 cube (c) vortex state 
in 1 x 1 x 1 state (c) homogeneously magnetized sphere with radius 0.5. 
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Method N relerr e relerr h av. relerr h [°] max. err h [°] 



DM 


40 x 40 x 40 


2.9e - 


- 09 














SP 


40 x 40 x 40 


Lie - 


-03 


Lie - 


-03 


2.3e - 


-05 


5.0e + 


00 


TG (CPr = 1) 


40 x 40 x 40 


3.8e - 


-04 


2.3e - 


-03 


6.9e - 


-06 


2.5e + 


00 


FES 


7.2e + 04 


8.6e - 


-04 


2.2e - 


- 03 


3.2e - 


-05 


5.2e + 


00 



Table 3: Homogeneously magnetized unit cube, relative error in the energy, the average relative 
error in the field/field- angle (w.r.t. DM). 



5 Numerical Experiments 

5.1 Homogeneously Magnetized Cube 

As a first benchmark we take a homogeneously magnetized unit cube and compute the magne- 
tostatic energy for varying grid-size N, where the exact value is e^ = 1/6 [hqM^]. Tab. 3 shows 
for each of the described methods the relative errors in the energy w.r.t. the exact value and 
the relative error in the field computed by (26), as well as the angular deviation (error in the 
field-angle) to the field computed by the FD-Demag method, see (10). We take the relative I2— 
error as a measurement for the field-error, i.e. 



relerr 



p=x,y,z 



\H 



(p) 

demag 



H 



(p) 

method 



1/2 



(26) 



The errors in the field- angle in Tab. 3 mostly occure at the edges of the cube. 



Fig. 7a shows magnetostatic energy calculations for different spatial discretizations. The 
DM method is almost exact and does not depend on spatial discretization. The reason is 
that the discretized demagnetization tensor is computed assuming homogeneously magnetized 
computational cells. Also the resulting stray field is analytically averaged per cell. Since the 
energy calculation is bilinear in the magnetization M and the stray field H, the error is a pure 
rounding error. 

The FES method is the slowest converging method for this problem. A possible reason is the 
large external stray field of this setup. The numerical integration of the diverging metric tensor 
g leads to an underestimation of the external space and consequently to an underestimation of 
the magnetic potential in the sample. Thus the FES method is particularly sensitive to setups 
with large external stray fields. 

SP and TG methods are also based on the computation of the scalar potential, whereby 
the field is obtained by finite differences. In [13] it is shown that the TG method essentially 
computes the scalar potential exactly for piecewise constant magnetization. The error in the 
energy in both methods (SP and TG) is mostly caused by numerical approximation of the gra- 
dient in the field computation, whereas TG shows the better approximation properties for this 
problem. In addition to it, TG uses an exact rank— 1 representation for the uniform magnetiza- 
tion which makes the computation sub-linear (namely 0(N 2 ^ 3 )) with small scaling factor and 
allows computations for dozens of millions cells without any problems related to storage and 
computational cost. 
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Method N e relerr h av. relerr h [°] max. err h [°] 



DM 


40 x 40 x 40 


1.528e 


-01 












SP 


40 x 40 x 40 


1.526e 


-01 


1.8e- 


03 


5.0e - 


05 


7.2e + 00 


TG (CP,r = 6) 


40 x 40 x 40 


1.529e 


-01 


1.8e- 


03 


7.8e- 


06 


2.6e + 00 


FES 


7.2e + 04 


1.526e 


-01 


2.5e- 


03 


6.0e- 


05 


6.8e + 00 



Table 4: Flower state for magnetization in the unit cube, energy, the average relative error in 
the field/field-angle (w.r.t. DM). 



5.2 Flower and Vortex State in a Cube 

We do the same comparison as in Sec. 5.1 for the flower state, see (27) and Fig. 5b, and Tab. 4 
for the results. The main magnetization direction is taken to be along the z- axis, and the flower 
is obtained through an in-plane perturbation along the y - axis and an out-of-plane perturbation 
along the x-axis. Assuming polynomial expressions for the perturbations, as in [8], our flower 
is the normalized version of 

™> x {r) = \xz, 

m y {r) = \yz + ^y 3 z 3 , (27) 
m z (r) = 1, 

where the center of the cube is located at (0, 0, 0). We choose a = c = 1 and 6 = 2. 

The results are similar to those of the homogeneously magnetized sample. In contrast the results 
of the DM method are not exact in this case, but Fig. 7b shows that the DM method converges 
faster than all other methods. 

The next comparison is for a vortex state in a unit cube, see Fig. 5c, described by the model 
in [14], i.e. 

m x {r) = - - (l-exp(-4-2)) 5 , 

2 

m y (r)= ^(l_exp(-4^))5, (28) 
m z (r) = exp ( - 2 -^), 

where r = x 2 + y 2 , and we choose the radius of the vortex core as r c = 0.14. The vortex 
center coincides with the center of the cube, and the magnetization is assumed to be rotationally 
symmetric about the x/y-axis and translationally invariant along the z-axis. The results can 
be found in Tab. 5. 

The most notable difference to the previous tests is the large field error in the FES method. 
It shows that the error occurs in the center of the vortex, where the gradient of the magnetization 
peaks. A possible solution for this problem would be an adaptive meshing, which is currenty 
not implemented in our FES code. 

5.3 Thin Film 

We first take a homogeneously magnetized 1 x 1 x 0.1 thin film (magnetization out of plane), 
see Tab. 6 for the results. Tab. 7 shows the results for the vortex state (out of plane) in the 
same thin film geometry. 
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Method N e relerr h av. relerr h [°] max. err h [°] 



DM 


40 x 40 x 40 2.189e- 


02 










SP 


40 x 40 x 40 2.163e- 


02 2.4e - 


- 03 


3.4e- 


04 


1.2e + 01 


TG (Tucker, 


r = 10) 40 x 40 x 40 2.193e- 


02 4.0e - 


- 03 


3.2e- 


04 


1.16 T" Ul 


FES 


7.2e + 04 2.160e- 


02 2.1e- 


- 02 


6.1e- 


02 




Table 5: Vortex state for magnetization in the unit cube, 
the field/field-angle (w.r.t. DM). 


ener 


gy, the averag 


;e relative error in 


Method 


N e 


relerr h 


av. relerr h 


:°] 


max. err /i [°] 


DM 


80 x 80 x 8 4.025e - 02 












SP 


80 x 80 x 8 4.021e - 02 


1.7e - 03 




2.6e - 05 




4.5e + 00 


TG (CP, r 


= 1) 80 x 80 x 8 4.025e - 02 


3.7e - 03 




6.4e - 06 




2.1e + 00 


FES 


4.9e + 04 3.983e - 02 


5.5e - 03 




1.9e-05 




5.0e + 00 



Table 6: Homogeneously magnetized 1 x 1 x 0.1 thin film (magnetization out of plane), energy, 
the average relative error in the field/field- angle (w.r.t. DM). 



The results for methods that do not rely on spatial discretization outside the sample perform 
equally well on this geometry. FES, instead, shows a deterioration of performance due to the 
worse ratio of shell and sample elements while leaving the number of DoF unchanged. 

5.4 Sphere 

As the last test a homogeneously magnetized sphere with radius R = 0.5 is simulated. The 
spatial discretization in case of cuboid grids is done by setting the magnetization M = (0, 0, M z ) 
in cells whose center lies within the sphere. This leads to staircase artifacts as shown in Fig. 8a. 
For the FES method the sphere is discretized such that the volume of the discretized sphere 
matches the anaytical volume. 

The magnetostatic energy for different spatial discretizations is displayed in Fig. 7d. The 
FES method shows the fastest convergence, which is obviously a consequence of the better 
approximation of the curved surface, see Fig. 8b. Also the field computation benefits from this 
better approximation, see Fig. 8c and 8d. 

Work was done on the treatment of curved surfaces within cartesian grid methods [9]. 
Still the use of irregular grids is a more natural way of describing curved surfaces and is thus 
preferable. 



Method 


N 


e 




relerr h 


av. relerr h [°] 


max. err h [°] 


DM 


80 x 80 x 8 


1.569e 


-03 








SP 


80 x 80 x 8 


1.555e 


-03 


2.4e - 03 


4.6e - 05 


3.6e + 00 


TG (CP,r = 1) 


80 x 80 x 8 


1.569e 


-03 


2.9e - 03 


1.8e-05 


4.0e + 00 


FES 


4.9e + 04 


1.496e 


-03 


6.0e - 03 


6.4e - 04 


2.1e + 01 



Table 7: Vortex magnetization in 1 x 1 x 0.1 thin film (magnetization out of plane), energy, the 
average relative error in the field/field-angle (w.r.t. DM). 
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Figure 8: Spatial discretization and stray field of a homogeneously, in z-direction magnetized 
sphere in the middle xz-plane (a) Finite difference approximation of the spherical sample with 
50 x 50 x 50 cells (b) Finite element approximation of the spherical sample with 9429 tetrahedra 
(including the shell elements) (c) z component of the stray field, calculated with DM method 
(d) z component of the stray field, calcualted with FES method 
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6 Conclusion 



We investigated several test magnetization configurations with different methods for the stray- 
field computation and compared the results. There is no clear winner in this comparison of 
numerical methods for the stray-field calculation. Computations on cuboid structures are best 
done with methods that compute on cuboid grids, namely the DM, SP and TG methods. The 
TG method is not only the fastest choice, it is also able to handle very large grids due to low- 
rank tensor approximation or representation of the magnetization. However the TG method is 
not yet well investigated in the context of full micromagnetic simulations. In order to preserve 
the sublinear complexity and storage requirement features further research on the behaviour of 
low-rank magnetization during energy minimization or LLG integration has to be done. 

The SP method is faster than the DM method by a factor of 1.5 and needs about 30% less 
memory. This speedup comes at the expense of accuracy. Among the Cartesian grid methods, 
the DM method is most accurate since the stray field is computed directly. Both the SP and 
the TG method show an additional error due to the finite-difference gradient computation. 

For curved structures FES is a good choice. The obvious reason for this is the use of irregular 
meshes, which are able to model the curvature much better than cuboid grids. In contrast to 
FEM/BEM methods FES can be implemented using sparse matrices only, since the presence of 
the dense boundary element matrix is overcome with the shell transformation. 
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A Low-Rank Tensor Formats 

A tensor A G R^x^xATg ig gaid tQ bg in canonica i f ormat (CANDECOMP/PARAFAC (CP) 
decomposition) with tensor product rank r, if 

r 

A = Aj uf ] o uf ] o uf ] (29) 
i=i 

with A; G R, vectors ttf G R^, and o is the vector outer product. A particular entry of a 
canonical tensor is given by 

a^^Xt^Uu^^ufX (30) 
i=i 

Abbreviating notation as in [22], a tensor A G rN 1 xN 2 xN 3 in cp format 

can be written as 

A = [A; U^\U^ 2 \U^}, (31) 

with weight vector A = [Ai, . . . , X r ] G R r and factor matrices = [u^ | • • • | ] G R N i xr . 
From (31) it can be seen that the number of degrees of freedom (DoF) of a CP tensor is 
r + r J2j Nj (compare with FJ ■ Nj for a dense tensor), also see Tab. 1. 
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A tensor A G R NlxN ' 2xN 3 [ s sa id to be in Tucker format (Tucker tensor) if it can be 
represented in the form 

A = C xil/i x 2 U 2 X3C/3, (32) 

with the so-called core tensor C G [R r i xr 2xr 3 anc j j ac ^ or matrices Uj G R-"i Xr i. 
Hereby the key-operation is the n-mode (matrix) multiplication of a tensor .4. G R Ar i xAr 2xAf 3 
with a matrix L7 G R xNn , which is the multiplication of each mode-n fiber of A by the matrix 
U, i.e. 

A x n U G R x ?=^, M,- = { ^' j * n (33) 

J I M, j = n. 

In contrast to CP tensors, the ranks in the Tucker representation can be different in each 
mode (dimension). In the discussions of Sec. 4 and the experiments in Sec. 5 we used the same 
rank r for each mode, i.e. r = rj. 

For a tensor in Tucker format Ylj r j + Ylj r j Nj entries have to be stored, which is a compression 
for rj <C Nj, also see Tab. 1. 

For a sum of Tucker tensors one can only store the factor matrices and core tensors of the 
summands, which is called block- CP format. 

Linear algebra operations for low-rank tensors, like the inner product, tensor-matrix product 
etc., can be performed without forming the dense tensors [3], which makes this operations faster 
than their conventional counterparts. 
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